Monte Carlo simulations with a generalized 
detailed balance using the quantum-classical 

isomorphism 

Yefim I. Leifman 
Department of Mathematics and Statistics, 
Bar-Ilan University, 
52900, Ramat-Gan, Israel 
leifmany@lycos.com 

February 2, 2008 

Abstract 

The main idea of this work is that the quantum-classical isomorphism is 
a suitable framework for a generalization of the notion of detailed balance. 
The quantum-classical isomorphism is used in order to develop a Monte Carlo 
simulation with controlled deviation from detailed balance, that is with a 
generalized detailed balance and known relative entropy with respect to the 
reference process at each point. In order to apply this method to molecular 
simulations a new algorithm for realization of a partial chirotope, based on 
linear programming, a new distance geometry algorithm and a new all-atom 
off-lattice Monte Carlo method are proposed. 
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1 Background 

The characteristic time of an event in the molecular world is 10 -15 sec, i.e. one 
iteration of a molecular dynamics must simulate changes corresponding to a time of 
this order. Biomolecules of interest, such as proteins, have thousands of atoms and, 
even using the simplest approximation for the molecular potential and a powerful 
computer, only thousand to million of iterations can be performed in a reasonable 
time. The approximate time of protein folding even in vivo is 10~ 5 sec. The timing 
for a standard Amber benchmark 159 residue protein in water is 249 ps/day of 
simulations on a single 3.4 GHz processor [2]. The performance of the NAMD 
program on different platforms can be viewed in [54] . However, significant changes of 
molecular conformations were achieved, for example, in all-atom molecular dynamics 
simulations of 36 residue protein on supercomputer of hundreds processors [TT] . 
Another method of molecular simulation is molecular Monte Carlo simulation. It 
does not try to simulate the physical movement of a molecule but only visit (sample) 
its conformational space according to an appropriate probability distribution, such 
as Boltzmann distribution. It gives a hope of closing the aforementioned time gap. 

A Monte Carlo (MC) simulation is a random process. Usually, it is a Markov 
process [19J. The Markov process {X t } igT is specified by its transition probability 
P(s,x,t,Y) = P(X t G Y\X S = x), s < t, and initial distribution - the distribution 
of the random variable X$. Y is in the smallest a-algebra that contains all open 
sets of the state space S of the process. If the transition probability depends only 
on the difference between s and t, that is, if there exists a function P(t, x, Y), such 
that P(s,x, t,Y) = P(t — s,x,Y), then the Markov process is called temporally 
homogeneous. A jump process is a continuous time process which changes its state 
after non-zero time. A temporally homogeneous Markov jump process is determined 
by its initial distribution, jump rate j(x) for a jump from each x G S and (jump) 
transition probability P(x,Y) with (jump) transition probability density p(x,y). 
Only temporally homogeneous Markov jump processes with j(x) = 1 for all x G S 
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will be considered in this work. In this case jumps occur according to the Poisson 
process. 

A process with the transition probability density p(x, y) satisfies the semi- detailed 
balance condition [13] if there exists a density /i such that for all y G S 



One says that such process with the initial density /x is in a steady state. 

A process with the transition probability density p(x, y) satisfies the detailed 
balance condition if there exists a density \x such that for all x,y G S 



One says that such process with the initial density \i is in equilibrium. 

If a process with a finite number of states satisfies the detailed balance condition, 
then the process converges to a limiting distribution, which is unique [26], [30]. The 
simplest example of a process with detailed balance is a process with independent 
outcomes, that is p(x, y) = p(z, y) and fi{x) = p(z, x) for all x,y,z G S. 

A Metropolis Monte Carlo simulation [2B] is an example of a process with detailed 
balance. Let given a system of particles in 3-dimensional space. Let r and n be 
vectors of coordinates of two states of the particle system in the phase space. The 
Boltzmann law gives the ratio of densities to be in the state r or n in equilibrium 



where E(r) is the energy of the conformation of the particles, T is a temperature, k 
is the Boltzmann constant. The aim of the Metropolis MC sampling algorithm is to 
sample the phase space according to the distribution which satisfies the Boltzmann 
law. First, one generates a conformation r (present). Next, one generates a new 
conformation n by adding a small random displacement to r. One must now decide 
whether the new conformation will be accepted or rejected. One wants to choose a 




/i(x)p(x,y) = fi(y)p(y,x). 



H{r) 



= exp 
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transition probability density p(r, n) such that the detailed balance condition 

n(r)p(r,n) = fi(n)p(n,r) 

is satisfied. We denote the probability density to try a move from r to n by l(r,n) 
and the probability of accepting a move from r to n by a(r,n). Assume that the 
transition probability density is given by 



One of the possibilities to satisfy this condition is the choice of Metropolis et al. [28] 




If E(n) < E(r), then the move is accepted. If E(n) > E(r), then we generate a 
random number U from the uniform distribution in the interval [0, 1] and we accept 
the move if U < a(r,n) and reject otherwise. I is not specified, except for the 
assumption that it is symmetric. This reflects freedom in the choice of moves. 

It is obvious that an infinite number of states of the gas corresponds to a given 
macroscopic condition of the gas. Through macroscopic measurements one should 
not be able to distinguish between two gases existing in different states (thus cor- 
responding to two distinct representative points in phase space) but satisfying the 
same macroscopic conditions. Thus when one speaks of a gas under certain macro- 
scopic conditions, one is in fact referring not to a single state, but to an infinite 
number of states. In other words, one refers not to a single system, but to a col- 
lection of systems, identical in composition and macroscopic condition but existing 



p(r,n) = a(r,n)l(r,n). 



If I is symmetric, i.e. l(r,n) = l(n,r) , then detailed balance implies 



n(r)a(r, n) = fi(n)a(n,r), 



and therefore, 
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in different states. Such a collection of systems is called an ensemble (Chapter 4 of 
|18j). which is geometrically represented by a distribution of representative points in 
phase space, usually a continuous distribution. An ensemble is completely specified 
by this distribution. Metropolis MC samples an ensemble which is called canonical 
ensemble. It is appropriate to a system whose temperature is determined through 
contact with a heat reservoir (Chapter 7 of [T8]). 



In generalized-ensemble simulations [33J, each state is weighted by a non-Boltzmann 




probability weight factor. This allows the simulation to escape from energy barriers 
and to sample much wider space than by conventional methods. Monitoring the 
energy in a single simulation run in such ensembles, one can obtain also canonical 
ensemble averages as functions of temperature. One of the best-known generalized- 
ensemble methods is the replica- exchange MC [23], [SS]- The system for a replica- 
exchange MC consists of non-interacting copies, or replicas, of the original system in 
canonical ensemble at different temperatures. There is a one-to-one correspondence 
between replicas and temperatures. A simulation of replica-exchange MC is then 
realized by alternately performing the following two steps. 

1. each replica in the canonical ensemble at a fixed temperature is simulated 
simultaneously and independently for a certain number of Metropolis MC 

steps, 

2. A random pair of replicas which are at neighboring temperatures T m and 
T m+ i are exchanged moving the low-temperature conformation to the high- 
temperature simulation and vice versa. These replica swaps are accepted ac- 
cording to the Metropolis criterion with the acceptance probability 



The effect of replica exchange is to prevent low-temperature simulations from be- 
coming trapped in local minima, because they are occasionally swapped to higher 
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temperatures where they can escape these minima and move to other regions of 
phase space. Simultaneously, the low-temperature simulations are always being 
seeded with low-energy conformations produced by simulations at higher tempera- 
tures. The replica-exchange MC sampling is in detailed balance. 

The notion of the quantum-classical isomorphism (sometimes it is called just the 
classical isomorphism) originated in Chapter 10 of [12] . The derivation can be found 
in Chapter 10 of pQ. Consider a quantum particle at temperature T, (3 = l/{kT). 
The density operator is p = exp(— (3H), where H = — + V{x). The density 
matrix is 

p(r,r',/3) =< r\exp(-/3H)\r' > 

=<r\ exp(-(3H/K) . . . exp(-(3H/K) . . . exp(-[3H/ K)\r' > . 

Inserting unity in the form 1 = j \r >< r\dr of the integral of projectors |r >< r| 
over the volume of the system between each exponential gives 

p(r,r',P) = J < r\ exp(-f3H / K)\r 2 >< r 2 \ exp(—{3H/K)\r 3 > 

... < r K -x\ exp(-(3H/K)\r K >< r K \exp(-f3H/K)\r' > dr 2 . . . dr K 
= J p(r, r 2 ,P/ K)p(r 2 , r 3 , 0/K) . . . p(r K , r' , p / K)dr 2 ■ ■ ■ dr K . 
If (3/K is sufficiently small, the following approximation is valid [12] 

p(r a ,r b ,(3/K) « p free (r a ,r b ,P/K) exp(-^(V d (r a ) + V d (r b ))) 

where V cl {r) is the classical potential energy, and the free-particle density matrix 
for a single particle of mass m is [12] 

p free (r a ,n^/K) = 2 exp(-^|d), 

where r ab = \r a — r&| and D is a dimension of the space. Therefore 

' (ri ' ri '#*(^) 2 / ex p(-i|^ + ^ + ---^i)) 



6 



exp 



(-jL(y d ( ri ) + V c \r 2 ) + ... V d {r K )))dr 2 ...dr K 



(1) 



The quantity p(r,r,j3) is proportional to the density to be in r. If an iV-particle 
system is considered, replace D by ND. That is, the density of the iV-particle 
quantum system which satisfies the Boltzmann law corresponds to the density of 
the .fTiV-particle classical system which satisfies the Boltzmann law and consists of 
K copies of the iV-particle classical systems and in the copies with neighbor numbers 
i and i + 1 for 1 < i < K and K and 1 the corresponding particles are connected 
by springs with spring constant mK/(/3 2 h 2 ). This approximation becomes exact as 
K — > oo [1] and can be used in a conventional MC simulation to investigate quantum 
properties. 

The main difficulties in all-atom detailed balanced MC simulations of biochem- 
ical processes involving big molecules are that these methods have huge autocorre- 
lation time (Section 2 of |40j) and these processes, generally, do not approach to 
equilibrium. 

A process is reversible if p(x, y) > implies p(y, x) > for all x,y 6 S. Let us 
mention some characteristics of non-equilibrium reversible processes. Such process 
satisfies the master equation 



dSa{^t)/dt is deduced from master equation and according to Section 2.4 of |23j 
it splits into the entropy production rate R(fit) and the entropy flow rate A(p t ) as 
follows 




For the probability density \it the Gibbs entropy is given by 




dS G (n t ) 
dt 



R(Ht) - A(fMt) 



where 



R (i*t) = a / Mx)p(x,y) - nt(y)p(y,x)) log 




x,y£S 



fit(x)p(x,y) 

th(y)p(v,x) 



dxdy 
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and 

A(jk)= [ {JH(x)I(x)dx = (I) Mt , I(x) = I p(x,y) \og P ^ V \ dy. 
Jx£S J yes P{y, x ) 

The entropy production rate is expressed in Section 3.1 of [15] in terms of the 
"particle fluxes" 

J t (x) = fit(x)p(x, y) - nt(y)p(y, x) 

and "forces" 

77? / \ , Ht{x)p{x,y) 

F t (x) = log—— -. 

fi t {y)p{y,x) 

As was mentioned in [6], entropy production rate can be considered as a measure of 
a lack of equilibrium. 

Definition 1.1 ( [1TJ ) For two probability distributions P and Q with probability 
densities p and q, the relative entropy (Kullback-Leibler divergence) is defined by 

D K l(P\\Q)= [ p(x)\og^\dx. 

If entropy is measured in bits, the logarithm in this formula is taken to base 2, or to 
base e, if entropy is measured in nats. 

The Gibbs inequality says that Dkl{P\\Q) > and the relative entropy is zero iff 
P = Q. The entropy flow rate A(fi t ) for the distribution fi t which is concentrated at 
one point x is the relative entropy of the transition probability density p'(y) = p(x, y) 
and the probability density p"(y) = p(y, x). 

2 The problems addressed in this work 

In this work we consider three different, but related problems: 

1. generalizing the detailed balance condition in order to include processes which 
do not preserve any distribution, but with the property that, roughly speaking, 
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it is known how much information the "jumper" must retrieve from its path in 
order to reach its current position (when the corresponding reference process 
is given), then building processes with such generalized detailed balance, 

2. enhancing all- atom off-lattice molecular MC simulations which are in detailed 
balance, for example replica-exchange MC, in order that they can be performed 
with a move set consisting of separate moves of each atom, with all degrees 
of freedom also in the case of dense atom packing; we do this by means of a 
new distance geometry algorithm, which plays in such MC simulations a role 
which is similar to the role of the SHAKE algorithm [36], [JT] in molecular 
dynamics, 

3. building an initial sample for all-atom off-lattice molecular MC simulations 
according to chirality constraints and distance constraints, and additionally, 
geometric manipulations with a molecule, which preserve, as far as possible, 
the aforementioned constraints, but exploit flexibility of a molecule. 

The connecting link of the following considerations of these problems is the 
distance geometry procedure which we call "centering". We build an initial sample 
of a molecule satisfying molecular chirality constraints and distance constraints with 
the help of linear programming and subsequent " iterative vibrant centering" . We 
perform a Metropolis MC with a move set consisting of separate moves of each atom 
in the sample space which is restricted with the help of " centering" . We combine 
this distance geometry method and this MC to exploit flexibility of a molecule in 
geometric manipulations with it. We use the same restricted sample space in MC 
without detailed balance. 

Section [3] provides an example of a process with generalized detailed balance 
with respect to its reference process using the quantum-classical isomorphism. 

Section H] provides the results of the numerical experiments considering some 
properties of the process which are similar to the example from Section [3J but can 
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be numerically examined. 

Clearly, there is a temperature at which the Amber force field cannot ensure 
the integrity of the molecule in the all-atom Metropolis MC sampling with a move 
set consisting of separate moves of each atom as described in [28]. This limits 
replica temperatures in replica-exchange MC. The same phenomenon can happen 
in non-equilibrium simulations. A way to overcome this difficulty is to change the 
potential in order to restrict such deformations. We change the potential with the 
help of distance geometry "centering" procedure. It can be considered as building 
the restricted sample space which includes all relevant conformations. It is described 
in Section [51 

In order to start a Metropolis MC with this potential one has to build an initial 
sample which is in the restricted sample space. As will be proved in Section [HI if one 
starts at some point of the sample space and performs as many centerings as needed 
from an infinite sequence of centerings which contains an infinite number of center- 
ings of each atom, then at certain step one reaches a point in the aforementioned 
restricted sample space. We call this algorithm "iterative centering". 

Molecular chirality constraints impose limitations on molecular conformations. 
These limitations are in addition to the limitations imposed by the weighted graph 
of the desired distances. A new algorithm for realization of a partial chirotope, 
based on linear programming is proposed in Section [7J 

Suppose, that a sample which satisfies a given partial chirotope (chirality con- 
straints) is built. Now we need to push this sample into the restricted sample space. 
Numerical tests show, that reiteration of iterative centering algorithm with chirality 
checking can become jammed if the initial sample is far from the restricted sam- 
ple space. For overcoming this difficulty the "vibrant iterative centering" distance 
geometry algorithm is proposed in Section [HI The vibrant iterative centering algo- 
rithm can be incorporated in a convenient computation scheme with a Metropolis 
MC in the restricted sample space. This scheme allows flexible manipulations with 
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a molecule with further equilibration. It is described in Section El 

Section [9] describes a molecular MC simulation without detailed balance using 
the quantum-classical isomorphism. In this simulation we use a process which is 
similar to the example which is described in Section EJ The simulation of Section 
El is only preliminary since the choice of appropriate parameters for true molecular 
simulations is not considered in this work, but the observations of Section [9] hints 
at the possibility to use such method in molecular MC. 

All algorithms with centerings are new. The iterative vibrant centering can be 
useful in existing distance geometry software in order to improve its sampling prop- 
erties. The restricted sample space can be useful in all-atom off-lattice molecular 
simulation software, for example, in order to increase the temperature of the hottest 
replica in replica-exchange MC. The algorithm for realization of a partial chiro- 
tope using linear programming is also new. It can be useful with iterative vibrant 
centering. 

The notion of generalized detailed balance in a framework of Langevin dynamics 
was proposed in [21]. As far as we know, our work is the first attempt to generalize 
the detailed balance condition in order to include processes which are not in a steady 
state but with the property, that it is known how much information the jumper must 
retrieve from its path in order to reach its current position (when the correspond- 
ing reference process is given). In our opinion the appropriate generalization of 
the detailed balance condition is the most important problem in computer molec- 
ular simulations, since we believe that an all-atom detailed balanced simulation of 
working ribosome will never be possible on a digital computer. 

The concepts of vertex, particle, bead and atom represent the same object in 
different contexts of this work. The notations of Sections El El d El are independent 
of the notations of Sections [H [31 BJ El 
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3 Generalized detailed balance 



Consider the Young's double-slit experiment (§§26,27,32,33 of [25]). Let two parallel 

plane screens be separated by a distance I. Let there be iV > 1 rectilinear parallel 

slits of width s on the first screen. In the case of the Young's double-slit experiment 

N = 2. Let the correspondent slit borders of neighbor slits be separated by a 

distance d. Let there be a plane monochromatic light wave with length A which 

propagates perpendicularly to the screens and hits the second screen through the 

slits in the first screen. If A <C s < d <C I, then the intensity on the second screen is 

_ J sin 2 Ysin 2 (iVgF) 
^ ~ iV 2 F 2 sin 2 (gY) ' 

where y is a distance between a point on the second screen and perpendicular pro- 
jections of the slits on the second screen, 

7TS 

Y = Ti y ' 

q = d/s > 1 and Jo is the maximum intensity on the second screen. In this case the 
Fraunhofer diffraction takes place, that is I (lyi) / 1 (I1J2) does not depend on /. 
Since 

°° sin 2 (ax) n 

— dx = a - 

x A 2 

(3.821.9 from [IB]), without loss of generality let the intensity distribution for N — 1 
be 

h(x) 



sin 2 x 



TCX 2 

Since 

50 sin 2 (ax) cos 2 (bx) . an 
-dx = — 



o 



x 2 4 

for < a < b (3.828.11 from [16]), let the intensity distribution for iV = 2 be 

j. , . sin 2 x sin 2 (2gx) 
27nr 2 sin 2 (gx) 

Then 

h(x) log 2 ^rldx = 1 
12 



since 



/ 

Jo 



oo 



log cos 2 (ax) 



■m 



(-l) n (6- 2an) 



cos(bx)dx 



an + Tib log 2 + 7r 



n 



where m < 6/ (2a) < m + 1, m = 0, 1, 2, 3, . . . (1.158 from [32]). 

In this example 1 bit of the information through what slit the photon passed 
translates in 1 bit of relative entropy Dkl(Ii{x)\\I2(x)) for all q — d/s > 1. The 
process "without information" can be considered as a reference process. The formula 
for Fraunhofer diffraction is an approximation, but this property can be directly 
verified in experiments on various screens. We call this property by the balance of 
relative entropy. 

In the path integral formulation of quantum mechanics [T2] the contribution 
of a particular path in the total probability amplitude for a photon has a phase 
proportional to time to travel along this path. In the aforementioned example the 
contribution of a particular path changes along this path such that the total intensity 
obeys the property of balanced relative entropy. A time with such property "rotates" 
in order to hide the information about the past which was not retrieved in time. 

The rest of this section is devoted to giving an example of the process with 
jumps which obey the property of balanced relative entropy. If the reference process 
is in detailed balance, then the condition of the balance of relative entropy can be 
considered as a generalization of the detailed balance condition. 

In order that the quantum-classical isomorphism be valid, the Boltzmann law 
must be satisfied [T2]. Therefore MC with independent outcomes of one quantum 
particle using the quantum-classical isomorphism with condition, that the vector 
of the coordinates a of the first copy of the particle is in the zero point must be 
defined by the following recurrent Levy construction [24] for 1 < n < K 



where each rj n is a vector with independent standard normal distributed coordi- 
nates. The Levy construction samples intermediate time points of a Brownian mo- 



K-n 



1 
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tion, conditioned to arrive to a predetermined point after a predetermined time. A 
Brownian bridge is a Brownian motion conditioned to return to the initial point 
after a predetermined time. The aforementioned MC sample (ao, . . . , a^-i) is the 
Levy construction for a Brownian bridge for the time interval divided into K equal 
subintervals with a = 0. In the case of the quantum-classical isomorphism, this 
Brownian bridge time parameter is considered as "imaginary time" [12] in contrast 
with ordinary time where the jumps of MC take place. In this work, imaginary 
time is discrete and denoted by a subscript and ordinary time is continuous, but a 
number of jumps in jump processes is denoted by a superscript. 

Consider MC simulations of two distinguishable quantum particles. Firstly we 
define an auxiliary process {N t }. 

Let (ao, . . . , dK-i) an d {bo, ■ ■ ■ , bx-i) be the aforementioned Levy constructions 
with ao = and bo = 0. Shift them by a and b according to the distribution 
of iVo, that is (aj, • • • , a#_i) = (ao, • • • , clk-i) + (a, . . . , a) and (&§, • • • , = 
(bo, . . . , b K _i) + (b, . . . , b). This is an initial sample for {N t }. 

Suppose, that i samples of {N t } are built, that is i — 1 jumps of {N t } have 
already happened. Build new Levy constructions (a , . . . , a t K _ 1 ) and (b , . . . , b l K _ 1 ). 
Randomly choose two numbers n\ and n\, n\ < n\, from to K — 1 with equal 
probability for each pair. If 

ll<i-^ill<ll<|-^ll (3) 

shift these new Levy constructions (a , . . . , a l K _ 1 ) and (b , . . . , by a l ni — a l ni 

and b l n% — b n , correspondingly with probability a, otherwise, shift them by a l ni — a l ni 
and b l , —b\. If II a % , — 6% II > ||a l i — b l 4 II shift them by a% — a% and b l t —b\ with 
probability a, otherwise, shift them by a* t — a % t and 6% — b % t and so on. 

{N t } satisfies (JTJ if it is in a steady state. Denote {N t } with a = 1 by {M t } 
and denote {N t } with a = 1/2 by {L t }. The process {L t } is a well-known MC 
of two free particles using quantum-classical isomorphism. Jumps of {M t } satisfies 
the balanced relative entropy condition with respect to jumps of {Lt}. Now, given 
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a process with aforementioned property and a detailed balanced process used as a 
reference process, we know how much information a jumper must retrieve from its 
path in order to reach its current position. 



4 Numerical experiments 

In order to present the results of the numerical experiments, we define some processes 
which are similar to {N t }. All these processes consider K copies of two particles 
in R. The copies are connected by springs as described in the discussion of the 
quantum-classical isomorphism in Section [TJ The first sample for all the processes 
of this section is ((0, . . . , 0), (8, . . . , 5)), where 5 > is large relative to the lengths 
of the steps of the processes. 

Define {N%} for < j < K as follows. Suppose, that i samples of {iVj?} are 
built. Build new Levy constructions (oq, . . . , a^^) and (b l , . . . , b> l K _]) according 
to ([2]). Randomly choose a number n\ from to K — 1 with equal probability 
for each number. Let n\ = n\ + j (mod K). If \\a % ni — b % n% || < \\a l ni — b % ni \\ shift 
these Levy constructions (a l , . . . , a % K _^) and (6q, . . . , b l K _ 1 ) by a % ni —a L n i and b % nl — b % n i 
correspondingly with probability a (these jumps we call "forward jumps") otherwise 
shift them by a 1 % — a 1 x andfo 4 ,— b \ ("backward jumps"). If \\a l , — b l , II > ||a%— &%|| 
shift them by a l n , — a l ni and b l ni — b ni with probability a ("forward") otherwise shift 
them by a 1 , — a\ and b % , — b % , ("backward") and so on. 

Define {W}} for < j < K as follows. Suppose, that i samples of {Wf} are 
built. Randomly choose a number n\ from to K — 1 with equal probability for 
each number. Let n % 2 = n\ + j (mod K) and h q = mK/(/3 2 h 2 ). 
If ||a l j — b l i || < 1 1 a 1 i — i||, then with probability a 



2 



>r> 1 — 1 1 n * 



1 





+ b\ 



n\+l 



1 



a 



2 



+ 



2 



+ 
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and other coordinates unchanged, otherwise 

, +1 _ fl nj-l + a »l+l 1 ,<+! _ 6 n' 2 -l + & n|+l _J_ * 

2 + ^^' ^ 2 + ^2^ 

and other coordinates unchanged. 

If || a 1 i — b l i || > || a 1 j — 6* 1 1|, then with probability a 



a't 1 = 1 1 ; to, &T = — 2 1 ; ?7 

2 ^WK ^ 2 

and other coordinates unchanged, otherwise 

a ^-i + a U+i i Ki-i + Ki+i i 

, + 1 _ " 2 1 "2 + 1 _|_ 71+1 _ "l 1 n l + 1 _|_ / 

* ~ 2 + V " 2 + 

and other coordinates unchanged. 

The following numerical experiments show that in some sense {W 7 /} and {iV/} 
are similar. Consider K = 8,16,32,64,128, a = 1,2/3,7/12,13/24,25/48. Let 
m be atomic mass unit, T = 300K. Standard normal distributed random num- 
bers are obtained by the Box-Muller algorithm [TU] from the uniformly distributed 
pseudo-random numbers [27] . We compute the following quantities: the num- 
ber of performed jumps J , the number T of forward jumps % < J for which 
\W i — b l i || — \\a % i — b l i || and lla't 1 — b l ^ l \\ — Wa 1 ^ 1 — fo't 1 !! have different signs, the 
number of such backward jumps 1Z, the average of the coordinates A = ■h J2n=i a n +1 
after the last jump, the average B = j^Y^,=i{^n +1 ~~ the average number 
C — Yun=i where C n is the number of i < J such that \\a % n — b l n \\ — \\a l n , — b l n ,\\ 
and ||a^ +1 — &^ +1 || — H^tt" — ^/^ll have different signs, where n' = n + j (mod X), 

(■A+BJATyOr , (A+B)J 

The results for {iV/} are shown in Table[T]of Appendix, the results for {W 7 /} are 
shown in Table [2] of Appendix. In the case of {iVj?} we take J jC = 2. The lengths 
in these tables are in units ^?=^- The observations are as follows. 

1. For both {N?} and {W 7 /}, if we fix K and j, then A and B are proportional 
to a - 1/2. 
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2. For {W 7 /}, if we fix K , then jrz^i is approximately constant. It is true also 
when j is not a constant, but a random variable, for example, uniformly dis- 
tributed in some interval, as for {N t } (Not shown in the Appendix). 

3. If we fix j/K, the closeness of t^z^; 101 {^t} an d for {W/} justifies that 
J7/C jumps of {Wf} approximate two jumps of {N}}. 

5 The restricted sample space 

Evidently, there is a temperature at which the Amber force field cannot ensure 
the integrity of the considered molecule in the all-atom Metropolis MC sampling 
with the move set consisting of the separate moves of each atom as described in [28J . 
This limits replica temperatures in the replica-exchange MC. The same phenomenon 
can happen in non-equilibrium simulations. A way to overcome this difficulty is to 
change the potential in order to restrain such deformations. 

A finite undirected weighted graph G is a triple < V,E,W >, where V denotes 
the set of its vertices, E denotes the set of its edges, and W : E — > M + is a 
function which specifies a positive weight for each graph edge. In order to restrict the 
sample space as was mentioned in Section [21 one has to set a weighted graph which 
corresponds to the molecule. Generally, atoms are vertices of this graph, covalent 
bonds form a part of its edges, pairs of atoms which are bonded through two covalent 
bonds form another part of its edges, weights are desired distances. The weights 
of edges which connect two atoms which are bonded through two covalent bonds 
determine the bond angles. So, the weighted graph of methane has 10 undirected 
weighted edges. The weighted graph of amide plane has also C a — C a and H — O 
edges, since their distances are well defined in amide plane. If one knows additional 
distances between atoms (for example from Nuclear Magnetic Resonance data), one 
adds corresponding edges and weights too. 

Let / : V — > M fc be a conformation of G in the /c-dimensional Euclidean space M fc 
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(regardless of the weights). Let the coordinates of vertices a v = f(v) be all distinct 
and let h : E — > IR + be the spring constants of edges. The point 

1 ,/r W({u, v}) . A ... 

will be called the center of the vertex (atom) u. 

The corresponding algorithm for finding the center of the vertex u takes as 
its input the adjacency-list representation [8] of the finite weighted graph G =< 
V,E,W> and current coordinates of its vertices. The adjacency- list representation 
of the graph consists of the array Adj of |V| lists, one for each vertex in V. For each 
u G V the adjacency list contains all the vertices v G V such that there exists 

an edge {u,v} G E. The weight W({u,v}) of the edge {u,v} G E is stored with 
a vertex v in it's adjacency list. /i({w, f}) are stored like W({u,v}) and the vector 
A[u] of current coordinates of the vertex u is stored with u. If x is a pointer to an 
element of the list Adjf-u], then, according to pseudocode conventions [8], vertex[x] 
denotes a vertex which adjacent to u ( denote it by v ), W[x] and H[x] denote 
W({u, v}) and h({u, v}). 

C enter (u) 

1 t <- 

2 g «- 

3 (0,0,0) 

4 /<-0 

5 x <— head(Adj[u\) 

6 while x ^ NIL 

7 do f <— vertex [x] 

8 2 <- A[u) - A[v) 

9 r <— ||z|| 

10 if r > 

11 then y <-y + H[x]{A[v] + (W[a:]/r)«) 
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12 t <- t + H[x] 

13 else f <- f + H[x]W[x] 

14 q <— q + 

15 x <— next [a;] 

16 if q > and t > 

17 then z «- - 

18 r <— ||2r|| 

19 if r > 

20 then y <— y + qA[u] + (f/r)z 

21 y ^(tl{t + q))y 

22 if t > 

23 then return (l/t)y 

24 else return 



As will be proved in Section [HI if one starts at arbitrary point of a sample 
space and performs as many centerings as needed from a sequence of centerings 
which contains an infinite number of centerings of each vertex, then at some step 
one achieves a point of a sample space such that for each u G V it holds that 
\\ a u — Cu|| < S(u) for a given > 0. Denote the set of such points by T>(S). 
As will be viewed in Section [6] the Hooke potential ^ egB jj (H e ll — W(e)) 2 is large 
outside T>(S). We change a given potential, for example the Amber force field, by 
assuming it infinitely large outside of T>(S) . S(u) has to be not too small if one does 
not want to restrict a considerable part of molecular degrees of freedom. 

6 The iterative centering algorithm 

Distance geometry is a part of computational geometry which is devoted to the 
study of the existence or non-existence of an embedding satisfying the condition in 
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the following definition as well as methods for construction of such embedding. 

Definition 6.1 Let G =< V, E,W > be a finite undirected weighted graph, where 
V denotes the set of its vertices, E denotes the set of its edges, and W : E — > IR + 
is a function which specifies a positive weight for each graph edge. An embedding of 
G in the k- dimensional Euclidean space M fc is a function f : V — > ¥L k such that for 
each edge e = {v,w} G E one has \\fiy) — f{w)\\ = W(e). G is called k-embeddable 
iff such an embedding exists. 

The problem of /c-Embeddability of an integer-weighted undirected graph is NP- 
hard [37]. However there is a semi definite programming algorithm for the Euclidean 
distance matrix completion problem, i.e. determining whether there exists a number 
k for which a given undirected weighted graph is fc-embeddable |22j. If there exists 
an embedding according to Definition 16.11 it is evidently in T>(S). 

Theorem 6.2 (|4J) A complete graph is k-embeddable iff each of its complete sub- 
graphs with k + 3 vertices is k-embeddable. 

Another distance geometry problem is that of bounded k-Embeddability, namely 
whether for given bounds l,u : E — > M + there exists a weight W with /(e) < W(e) < 
u(e) for which a graph G =< V,E,W> is fc-embeddable. 

There are number of methods which can be applied to solving the aforementioned 
problems, which generally arise in Nuclear Magnetic Resonance data interpretation: 
metric matrix distance geometry [5], [31], simulated annealing, variable target func- 
tion optimization [7J and global continuation [29J. 

The iterative centering distance geometry algorithm consists of performing as 
many steps as needed from an infinite sequence of centerings which contains an 
infinite number of centerings of each vertex. 

Proposition 6.3 Let {uk} be an infinite sequence of vertices of G. Let us apply a 
sequence {C enter (uk)}i<k< n of n centerings. Denote the displacement of the vertex 
u n after the iteration Center(u n ) by l n . Then lim^oo ||/ n || = 0. 
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Proof. Without loss of generality suppose that all coordinates of u n after n — 1 
iterations are zero. Denote the coordinate vectors of vertices adjacent to u n after 
n — 1 iterations by a%, a m , the corresponding weights by W\, ...,w m and the spring 
constants by hi, h m > 0. If ||aj|| > for every i,l < i < m, then according to (j4]) 
the coordinates of u n after n iterations is 



y 



a,-. 



Therefore 



( E h ) M 2 = H h \{ 1 - SHI - E HI ( x - jS, H - f 



i=i j=i " "" i=i 

similar to the Huygens theorem about momenta. 
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by the triangle inequality and 
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Therefore 



m m 
i=l i=l 



rn 



i=i i=i 



EMIk - 2/11 - ^) 2 - 



(5) 



If ||aj|| > for every i, 1 < i < k, and ||aj|| = for every i, k + 1 < i < m, we 
take some z, \\z\\ = 1 and we define 

l/= ^m r ( EM 1 ~ irnr) a » + E hi 

e—dl— 1 j—^ ii n 2=fc+l 

Similarly to the previous, we have (jSJ). We put 

XX 1 



«i 2=1 



as in the Center(u) algorithm. The right hand side of (jSJ) is twice the difference of 
the old and the new values of the Hooke potential J2 e ^E ^'(IMI — W{e)) 2 . That 
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is, at stage n, the Hooke potential decreases by at least ( IKJ| 2 > where l n is 

the shift of the center. If \\l n \\ 0, then the Hooke potential would become negative 
at some n. Since Hooke potential is non-negative we have lim n ^ 00 = 0. □ 

7 The realization of the partial chirotope related 
to the molecular chirality constraints 

Let xi, X2, ■ ■ ■ ,x r , yi, |/2, . . . , y r £ K r - Then the following Grassmann-Plucker rela- 
tion holds 

det(xi, x 2 , . . . , x r ) ■ det(yi, y 2 , . . . , y r ) 

r 

= det(yi, x 2 , • • • , x r ) ■ det(j/i, . . . , x x , . . . , y r ). 

i=i 

The difference of the left and the right sides is an alternating multi-linear form in the 
r+1 arguments xi, yi, yi-, ■ . . , y r , which are vectors in an r-dimensional vector space; 
hence, the difference of the left and the right sides is identically zero. For example, in 
rank 3 one gets for every set of 5 vectors (denoting determinants by square brackets, 
and labeling the points 1 to 5) the relation [123] [145] - [124] [135] + [125] [134] = 0. 
This requires that these 6 signs of the brackets on the left side are such that the 
equality is at least possible for this sign pattern, when actual scalars are not given: 
for example, these 6 signs could be +,+,+, +,+,+, but not +, +,-,+,+, +. 

Definition 7.1 ([3]) Let r > 1 be an integer, and let E be a finite set. A chirotope 
of rank r on E is a mapping x '■ E r — ► { — 1,0,1} which satisfies the following 3 
properties: 

1. X is not identically 0, 

2. x is alternating, that is, x( a vi, a <r 2 > ■ ■ ■ > a °>) = s ig n ( cr )x( a i; a 2, ■ ■ ■ , a>r) f° r a ^ 
ai, a,2, ■ ■ ■ , a r G E and every permutation a, 
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3. for all ai,a 2 , ■ ■ ■ , a r , bi, b 2 , ■ ■ ■ ,b r e E such that 
X(ai, 02, . . . , a r ) ■ x(h, 6 2 , . . . , b r ) ^ 0, 
t/iere exists an % e {1, 2, . . . , r} suc/i £/ia£ 
X(6i, a 2 , . . . , a r ) • xOi, • • • , ai, 6»+i, • • • A) 
= X(ai,a 2 , . . . ,a r ) • x(^i^2, • • .,6 r ). 

The axioms comes from abstracting sign properties in the Grassmann-Plucker 
relations for r-order determinants. The chirotope axioms are a version of the oriented 
matroid axioms. 

Suppose E = {1,2, ... ,n}. Given any (n — r)-tuple ( ) of elements 

in E, then we write (a[, . . . , a' r ) for some permutation of E\{a\, . . . , a„_ r }. Then 
{ax, ... , a n _ r , a[, . . . , a' r ) is a permutation of (1,2, ... ,n), and we can compute 
sign(ai, . . . , a„_ r , a[, . . . , a' r ) as the parity of the number of inversions of this string. 
The mapping x* '■ E n ~ r — > { — 1, 0, 1}, defined by 

X*(ai, ■ ■ ■ , On-r) = xK, • • • , aj.)sign(ai, . . . , a„_ r , a[, . . . , a' r ) 

is called the chirotope dual to the chirotope x- 

Let x '■ E r — > { — 1,0,1} be a chirotope, E = {l,...,n}. If there exists 
{xx, . . . , x n } C M r such that 

X(ai,a 2 , ...,a r ) = sign(det(x ai , x a2 , . . . ,x ar )) 

for all 1 < ai < a 2 < . . . < a r < n, then x is called realizable and q : E — > W, 
i i— > is called a realization of x- 

Let G r (R ra ) be the real Grassmann manifold of r-dimensional linear subspaces in 
]R n , or equivalently Mat rxn (R)/GL r (R), which corresponds to the space of config- 
urations of n vectors in R r modulo the action of the general linear group GL r (M). 
Thus the realization q of x corresponds to a point in G r (M ri ). The set of such points 
is called the realization space of x- 
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Let x '■ E r — > {—1, 0, 1} be a chirotope. Then for each subset {a 1; . . . , a r+2 } of 
E there exists {xi, . . . , x r+ 2} C M r such that 

X(a h ,a i2 , . . .,a ir ) = sign(det(x il , x h , . . .,x ir )) 

for all 1 < i\ < %i < . . . < i r < r + 2. This feature of chirotopes is called by local 
realizability. Local realizability follows from the facts that realizability is preserved 
under duality since G r (W l ) = G n - r (M. n ) and that all rank 2 chirotopes are realizable. 

The realizability problem for chirotopes is NP-hard [3]. There is an algorithm 
for a realization of a chirotope x which is valid when the realization space of x is 
contractible and x '■ E r — > {—1, 1} [2], [9]. 

If the alternating map x is only partially defined and the Grassmann-Plucker 
relation holds whenever x is defined on all its participants, then x is called a partial 
chirotope. A partial chirotope x' °f rank r on E is called extendable if there exists a 
chirotope % °f rank ron£ and for any ai, . . . ,a r G E, x'( a ii ■ ■ ■ > a r) — x( a i? • • • , a r ) 
holds whenever x'( a i> . . . , a r ) is defined. The problem of testing extendability of a 
partial chirotope is NP-complete [4"2] . 

Molecular chirality constraints impose limitations on molecular conformations. 
These limitations are in addition to the limitations imposed by weighted graph of 
desired distances. A set of inequalities of type det(xb — x a ,x c — x a ,Xd — x a ) > 
0, where a, b,c,d £ V, corresponds to molecular chirality constraints. Then the 
corresponding equalities x{ a i b, c,d) — 1 define a rank 4 partial chirotope. If o ( — > 
(x\,Xa,xl) satisfies these inequalities, then a \— > (l,x\,x^,x^) is a realization of 
the corresponding partial chirotope. The most widely adopted method to realize 
a partial chirotope related to molecular chirality constraints is the minimization of 
the function, which includes deviations from given oriented volumes, by simulated 
annealing starting from an approximate embedding [31]. An example of a realization 
of a partial chirotope by use of such function can be found in |43j . 

Let "maximize cx with conditions Ax < b and x > 0" be a linear program, x > 
means Xj > for all j. Particularly, let Cj be the price per unit of the product j 
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produced, Xj be the quantity of the product j produced, 6j be the quantity of the 
material % on hand, be the quantity of the material % required to produce one unit 
of the product j. Let yi be the price per unit of the material i. One is interested in 
selling the materials instead of the products if A T y > c. The dual linear program 
"minimize by with conditions A T y > c and y > 0" answers the question what is the 
minimal price of all materials when it is advantageous to sell the materials instead 
of to work. This price is the same as the maximal income in the first (primal) linear 
program. It is the figurative formulation of the linear programming strong duality 
theorem as economists learn it. 

In some cases the following algorithm allows one to realize a given molecular 
partial chirotope. Let Y C V be a set of vertices, whose coordinates appear in 
inequalities of type det(xb — x a ,x c — x a ,Xd — x a ) > 0. Without loss of generality 
one can demand det(:rb — x a , x c — x a , Xd — x a ) > e > for all these inequalities and 
Zi = xf > for all i 6 Y . If we fix x\ and xf for all % 6 Y then the inequalities 
become linear. A feasible solution of the following (symmetric) linear programming 
problem 

Minimize z\ + ... + z m subject to Az > e and z > 
is a solution of our problem. Its dual problem 

Maximize et± + ... + et^ subject to A T t < 1 and t > 

has zero as a feasible solution. If the original problem has a feasible solution, then its 
dual is bounded by the strong duality theorem. The dual problem can be solved by 
the primal simplex procedure and if it is bounded, then the solution of the original 
problem can be taken from the last simplex tableau, according to the Chapter 4 of 

In this method one has to fix x\ and xf for all i 6 F. Since e is a minimal 
volume for a parallelepiped spanned by Xb — x a ,x c — x a ,Xd — x a for each ordered 
base (a,b,c,d), the points (x a ,x1), (x\,x\), (xl,xf), (x d ,x^) cannot be on the same 
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straight line. We place all atoms of a molecule in a sequence and choose x\ = 
cos(27T2/n), x\ = sin(27ri/n) for all i G V and e = (sin(27r/n)) 3 , where n is a number 
of atoms. 

In practical implementation of this algorithm of realization of a molecular partial 
chirotope, one has to set a partial chirotope of a given molecule. For example, for 
C a atom of amino acid residue it is necessary to set 3 ordered bases. Fixing only 
2 of them jams vibrant iterative centering algorithm (the modification of iterative 
centering which will be described in Section |H]) and fixing 4 of them is too restrictive 
for the choice of x\ and xj. The fourth ordered base will be recovered by means 
of distance constraints. Similarly, for C a atoms of one spire of a-helix, in which 
participate 5 residues, it is necessary to set 3 ordered bases. This partial chirotope 
will be used also for chirality checking CheckChirality(u). 

Consider an example of poly-L-threonine Thr 180 . Each Thr residue contains two 
chiral centers. For C 13 atom of Thr residue it is necessary to set 3 ordered bases. 
Arrange 14 atoms of each Thr residue in a following order H-N-H-C a -C^-H-0 7 -H- 
C 7 -H-H-H-C-0 (or in the notations of Protein Data Bank H-N-H-CA-CB-H-OG1- 
H-CG2-H-H-H-C-0). The proposed algorithm successfully finds a realization of a 
corresponding partial chirotope. Let us add to this chirotope also the constraints 
on C a atoms which appear assuming Thrigo is twisted in 50 spires of right a-helix. 
The algorithm successfully finds a realization in this case. 

8 Metropolis Monte Carlo in the restricted 
sample space 

If a given partial chirotope is realized, one has to transform this realization, keeping 
correct chiralities, in order to achieve T>(S) and then to start the Metropolis MC 
simulation in the restricted sample space. Let C heckC hirality(u) be a function 
which checks whether quadruples of vertices which contain a vertex u satisfy a given 
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partial chirotope. Numerical tests show that reiteration of the iterative centering 
algorithm with chirality checking can become jammed if the initial sample is far from 
T>(S). (For example, consider the weighted graph with four vertices and four edges 
on a plane: let the starting configuration be A — (0, 0), B — (4,3), C = (4,-3), 
D = (0,5), the weights of (A,B) and (A,C) be 5, the weight of (B,C) be 6, the 
weight of (A,D) be 0.01, B,C,D be counter-clockwise and S(u) = 0.001 for all 
vertices.) For overcoming this difficulty one can use the following modification of 
the iterative centering algorithm: 

Vibr antC enter (u) 

1 a <- A[u] 

2 z <— Center{u) 

3 r «— \\a — z\\ 

4 if r > C ■ S[u] 

5 then A[u] <— a + C ■ S[u] ■ ((z — a)/r + c - RandomVectorQ) 

6 else if r > S[u] 

7 then A[u] ^— a + S[u] ■ ((z — a)/r + c - Random]/ ~ector()) 

8 else A[u] <— z + S[u) ■ c ■ RandomVectori) 

9 if not CheckChirality{u) 

10 then A[u] <— a 

RandomVectorQ denotes a function, which returns a uniformly distributed vec- 
tor in a sphere of radius 1 with the center at the origin. A coefficient c > 1 is 
introduced for ergodicity. A coefficient C > 1 is introduced for speeding-up. There 
is no guarantee that T>(S) will be achieved, but the examples of this section show 
that the method works. 

The vibrant iterative centering algorithm can be useful if the method described 
in Section [7] fails. Split some vertex from Y into several vertices and spread the 
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inequalities of the form det(x b — x a , x c — x a , x^ — x a ) > e in which the coordinates of 
the initial vertex participate over these new vertices. In the weighted graph put the 
desired distances between the new vertices be 0. Do this for several vertices from Y . 
Apply the described linear programming method and the vibrant iterative centering 
to bring nearer the vertices obtained from the same vertex. 
Introduce CheckDistance function: 

CheckDistance(u) 

1 x <— head(Adj[u\) 

2 while x ^ NIL 

3 do v ^— vertex[x] 

4 if \\A[v] - Center(v))\\ < S[v] or S[v] = 

5 then x <— next[x] 

6 else return FALSE 

7 return TRUE 



Let us sum up the proposed methods in the following computing scheme: 



TrialMove{u) 

1 a <- A[u] 

2 A[u] a + S[u] ■ RandomVectorQ 

3 if \\A[u] — Center(u)\\ < S[u] and CheckDistance(u) and CheckChirality(u) 

4 then E <— potential(u) 

5 z <— A[u] 

6 A[u] <— a 

7 e <— potential(u) 

8 if E < e or RandomQ < exp((e - E)/(kT)) 

9 then A[u] <— z 



10 else A[u] <— a 



28 



11 z <— Center(u) 

12 r <— ||a — z|| 

13 if r > S^w] or not C heck Di stance (u) 

14 then if r > C ■ S[u] 

15 then A[u] <— a + C • £>[«] ■ ((z — a)/r + c - RandomVectorQ) 

16 else if r > S'ftt] 

17 then A[u] <— a + Sf-u] • ((z — a)/r + c ■ RandomVector()) 

18 else A[u] <— z + S[u] ■ c ■ RandomVector() 

19 if not CheckChirality(u) 

20 then A[u] <- a 



RandomQ generates a uniformly distributed in [0, 1] random number [27] . 

Now we apply the proposed methods in their natural succession: firstly a real- 
ization of partial chirotope, then vibrant iterative centering and then the Metropolis 
MC in the restricted sample space. If a conformation which satisfies chirality con- 
straints is not given, then one can use the method of Section [7] to build such confor- 
mation. It is useful to rescale this conformation to its natural scale and proportions. 
Then one can start vibrant iterative centering with large S[u] to break frozen parts 
of the initial conformation and gradually decrease S[u] to required values. Then it 
is possible to start TrialMove over all atoms, which drives a molecule to T>(S) and 
then becomes to be the Metropolis MC simulation in T>(S). Similar to the original 
Metropolis MC [58] . only coordinates of one atom are changed near its current po- 
sition in a trial move. It makes the all-atom (including hydrogen) Metropolis MC 
possible also in the case of dense atom packing. Also TrialMove allows flexible 
manipulations with molecule by adding or removing weighted edges of the weighted 
graph and subsequent equilibration. 

Consider an example of poly- L- alanine Ala36, which is twisted in 10 spires of 
right a-helix. We add to the weighted graph, which is derived from the primary 
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structure of the molecule, distances of hydrogen bonds 0^ — — N( i+4 ) 

and distances C?k — C? i+2 \, which are well defined in a-helix (the parentheses contain 
numbers of residues). Then we apply the simplex procedure, vibrant iterative cen- 
tering algorithm and Metropolis MC simulation in the restricted sample space using 
Amber force field as described in the previous sections and receive the expected 
structure. 

9 Molecular MC simulation without detailed 
balance using the quantum-classical 
isomorphism 

In order to proceed to non-equilibrium molecular simulations we add pairs of parti- 
cles similar to that described in the example of Section [3] to the considered molecule. 
These pairs of particles are used as artificial devices and do not represent physical 
particles. In this section we shall call these artificial added particles by beads for 
convenience. Also we add one Hooke term per bead to the molecule potential so 
that it connects a bead to some atom of a molecule by a spring with spring constant 
h a and zero length when the spring is relaxed. Suppose, that a sample which satis- 
fies a given partial chirotope and distance constraints of such equipped molecule is 
built and equilibrated by methods described in previous sections. Subsequently we 
produce K copies of this system which are connected by springs as described in the 
discussion of the quantum-classical isomorphism in Section [TJ 

In this loaded case we cannot use two Levy constructions for an added pair of 
beads as described in Section [31 but in order to approach to a process whose jumps 
obey the property of the balance of relative entropy let K = 2j and choose the 
Z-th added pair of beads with probability pi, then uniformly choose two copies n\ 
and n\ = n\ + j, fix the n^-th copy of one bead from the /-th added pair and the 
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n l 2 -th copy of another bead from the Z-th added pair according to ([3]) and perform 
predefined large number of Metropolis MC steps over the rest of atoms and beads 
as described in Section [8] and so on. 

Consider the example of the linear polymer molecule which contains N = 16 
identical atoms with some Lennard- Jones constants and with neighbor atoms con- 
nected by springs. We add one aforementioned pair of beads to each pair of neighbor 
atoms. We constrain the sequence of the second beads of the added pairs to have 
chiralities of right helix and produce two copies of this system which are connected 
by springs as described in Section [TJ Then we start the simulation and observe that 
the polymer moves ahead. If one fixes the last atom in the space, then the polymer 
twists around the fixed atom like boa. The twisting polymer squeezes itself out. 
This observation hints on the possibility to use such method in molecular MC. 
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Appendix 



The following notations are used in Tabled] and Table [2] (see Section 0] for details). 
K the number of copies in the quantum-classical isomorphism, 
j the superscript parameter of {iVf } or {W/}, 
a the probability parameter of {iV/ } or {W 7 ^}, 
J the number of forward and backward jumps, 

JF the number of forward jumps % < J for which II a 1 1 — b % t II — II a 1 4 — b l 4 II 

"l n l n 2 n 2 

and lla't 1 - ft*! 1 II - lla't 1 - ft*! 1 II have different signs, 
7£ the number of backward jumps i < J ioi which II a 1 * — b % 4 II — ||a% — 6% II 
and H^t 1 — II — l|a l t — fr l t II have different signs, 

ii nj n| " 11 n 2 n 2 

A the average of coordinates A = J2^=i a n +1 after the last jump, 

B the average B = -j^2^^=i(^n +1 — 5), where the first sample for the 

processes is ((0, . . . , 0), (5, . . . , 5)), 
C the average number C — 4 J2n=i ^n, where C n is the number of i < J 

such that \\a l n - b l n \\ - \\a % n , - b % n ,\\ and ||a^ +1 - - Ha^t 1 - fr^T 1 !) have 

different signs, where n' = n + j (mod K), 
In the case of {iV/} we take J7"/C = 2. The lengths in these tables are in units 

h 

VmkT' 
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K 


j 


a 


J 


T 


n 


A 


B 


(A+B)J 
(T-IVjC 


8 


2 


1.0000 


10000000 


5002391 





1.2212c+006 


1.2242e+006 


9.778e-001 


8 


2 


0.6667 


10005000 


3336249 


1668173 


4.0719c+005 


4.0778c+005 


9.773c-001 


8 


2 


0.5833 


10000000 


2919143 


2081936 


2.0472c+005 


2.0419e+005 


9.768e-001 


8 


2 


0.5417 


20000000 


5415950 


4583516 


2.0225c+005 


2.0394e+005 


9.758e-001 


8 


2 


0.5208 


20000000 


5205748 


4792723 


l.Olllc+005 


9.7748c+004 


9.629c-001 


16 


1 


0.6667 


10000000 


3331214 


1667291 


2.2680c+005 


2.2837c+005 


5.471c-001 


16 


2 


0.6667 


10000000 


3333718 


1666593 


3.1050c+005 


3.1231c+005 


7.470c-001 


16 


4 


0.6667 


10000000 


3332375 


1668024 


4.0726e+005 


4.0487e+005 


9.763e-001 


16 


8 


0.6667 


10000000 


3332983 


1667959 


4.7057c+005 


4.6881e+005 


1.129c+000 


32 


1 


0.6667 


10000000 


3333992 


1665310 


1.6566c+005 


1.6342e+005 


3.946e-001 


32 


2 


0.6667 


10000000 


3334183 


1666927 


2.2823e+005 


2.2691c+005 


5.461c-001 


32 


4 


0.6667 


10000000 


3335290 


1664622 


3.1246c+005 


3.1090e+005 


7.460c-001 


32 


8 


0.6667 


10000000 


3334115 


1665497 


4.0766c+005 


4.0893c+005 


9.788e-001 


64 


2 


0.6667 


10000000 


3333389 


1665903 


1.6270e+005 


1.6609c+005 


3.944e-001 


64 


8 


0.6667 


10000000 


3334884 


1667409 


3.1174c+005 


3.1077c+005 


7.465c-001 


64 


16 


0.6667 


10000000 


3333361 


1665149 


4.0878e+005 


4.0681e+005 


9.783e-001 



Table 1: Numerical results for {iV/} 
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K 


j 




a 


J 


T 


71 




A 




B 


C 


(A+B)kVk 




16 


1 





.6667 


10000000 


2294496 


1153377 


1 


,5090e+004 


1 


5224c+004 


512272 


1 


700e+000 


5.187e-001 


16 


2 





.6667 


10000000 


1642125 


827546 


1 


,1095e+004 


1 


1140e+004 


369874 


1 


.747e+000 


7.381c-001 


16 


4 





.6667 


10000000 


1243440 


631437 


8 


4986e+003 


8 


5709e+003 


282097 


1 


.785e+000 


9.887e-001 


32 


1 





.6667 


10000000 


2251225 


1130260 


5 


,2767o+003 


5 


2874c+003 


251470 


1 


706e+000 


3.748e-001 


32 


2 





.6667 


10000000 


1580763 


798450 


3 


7948e+003 


3 


7869e+003 


178242 


1 


.754c+000 


5.436e-001 


32 


3 





.6667 


10000000 


1307517 


663010 


3 


1831e+003 


3 


1756e+003 


148053 


1 


.786e+000 


6.665e-001 


32 


4 





.6667 


10000000 


1138660 


582221 


2 


7568e+003 


2 


7539B+003 


129549 


1 


,792e+000 


7.644c-001 


32 


5 





.6667 


10000000 


1042861 


534477 


2 


,5251o+003 


2 


5549e+003 


118876 


1 


809e+000 


8.405c-001 


32 


6 





.6667 


10000000 


967434 


495375 


2 


,3921o+003 


2 


3670e+003 


110380 


1 


.825c+000 


9.131c-001 


32 


7 





.6667 


10000000 


911386 


467753 


2 


,2203e+003 


2 


2288e+003 


104017 


1 


.815c+000 


9.644e-001 


32 


8 





.6667 


10000000 


874623 


449650 


2 


,1260e+003 


2 


1291e+003 


100035 


1 


812e+000 


l.OOlc+000 


32 


g 





.6667 


12000000 


1024333 


527608 


2 


,5445e+003 


2 


5039e+003 


117076 


1 


,840e+000 


l.()41e+000 


32 


10 





.6667 


11000000 


911136 


470663 


2 


,2247e+003 


2 


2062e+003 


104318 


1 


.821e+000 


1.061e+000 


32 


n 





.6667 


14000000 


1142519 


589629 


2 


,7743o+003 


2 


7949c+003 


130899 


1 


.823c+000 


l.()77c+000 


32 


12 





.6667 


11000000 


874736 


453702 


2 


,1478o+003 


2 


1269e+003 


100289 


1 


838e+000 


1.114e+000 


32 


13 





.6667 


14000000 


1102885 


571297 


2 


7163e+003 


2 


7()51c+003 


126588 


1 


,846e+000 


1.128o+000 


32 


14 





.6667 


12000000 


938220 


487280 


2 


,2964e+003 


2 


2897e+003 


107856 


1 


.841c+000 


1.131c+000 


32 


15 





6667 


14000000 


1086256 


568979 


2 


,6827e+003 


2 


6788e+003 


125154 


1 


.876e+000 


1.159e+000 


32 


1 


1 


.0000 


10000000 


3332721 





1 


,5819e+004 


1 


5825e+004 


247677 


1 


719e+000 


3.834e-001 


32 


2 


1 


.0000 


10000000 


2369493 





1 


,1341o+004 


1 


1354e+004 


179076 


1 


.734e+000 


5.347c-001 


32 


2 





.5833 


10000000 


1381226 


991383 


1 


,9456o+003 


1 


8973c+003 


178013 


1 


.785e+000 


5.536e-001 


32 


2 





.5417 


20000000 


2573137 


2181102 


1 


,9085o+003 


1 


9266e+003 


356466 


1 


.771c+000 


5.491c-001 


32 


2 





.5208 


40000000 


4950753 


4559821 


1 


,9429o+003 


1 


8834c+003 


712677 


1 


.772e+000 


5.496e-001 


64 


1 





.6667 


10000000 


2233592 


1123461 


1 


8555e+003 


1 


8491c+003 


124756 


1 


.708e+000 


2.675e-001 


64 


2 





.6667 


20000000 


3102155 


1571597 


2 


6385e+003 


2 


6427e+003 


175391 


1 


.767e+000 


3.935e-001 


64 


4 





.6667 


20000000 


2208690 


1127635 


1 


,9124o+003 


1 


9001e+003 


125776 


1 


806e+000 


5.610e-001 


64 


8 





.6667 


22000000 


1782288 


917200 


1 


,5424e+003 


1 


5407c+003 


101851 


1 


.825c+000 


7.699e-001 


64 


16 





.6667 


30000000 


1909428 


988772 


1 


,6775o+003 


1 


6584e+003 


109624 


1 


.855c+000 


9.917e-001 


64 


24 





.6667 


25000000 


1450673 


755564 


1 


2598e+003 


1 


2479c+003 


83414 


1 


.847c+000 


1.081B+000 


128 


1 





.6667 


10000000 


2221748 


1119251 


6 


5420e+002 


6 


,5077e+002 


62092 


1 


714e+000 


1.906e-001 


128 


16 





.6667 


10000000 


573255 


297429 


1 


,7767e+002 


1 


7490e+002 


16466 


1 


.851e+000 


7.764e-001 



Table 2: Numerical results for {W(} 
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